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Abstract 

A  laboratory  sediment  sample  is  digitized  by  X-ray  computed  tomography  system  to  obtain  a  host  matrix  of  size  100^  of 
porosity  0.282.  Flow  through  the  porous  medium  is  studied  by  an  interacting  lattice  gas  model  for  fluid  driven  by  pressure 
bias.  The  bias  H(Q^H ^  1)  is  implemented  via  soft  and  hard  pressure  constraints  in  selecting  the  stochastic  movement  of 
fluid  particles.  The  response  of  fluid  flux  density  to  bias  H  is  nearly  linear  with  soft  pressure;  it  becomes  nonmonotonic 
with  the  hard  pressure 
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Understanding  the  flow  through  porous  media  has  attracted  a  considerable  interest  [1-4]  for  decades  due  to 
its  complex  response  properties.  Some  applications  include  secondary  oil  recovery,  spreading  of  spilled 
contaminants,  and  dissociation  of  gas  (methane,  carbon  monoxide,  etc.)  in  natural  occurring  processes  [5]. 
Heterogeneous  porous  media  for  geomarine  applications  are  generally  modeled  by  percolation  [6]  and  its 
variants.  Transport  and  flow  through  such  porous  matrices  have  been  extensively  studied  near  and  above  the 
percolation  threshold  [7,8].  Scaling  behavior  of  such  response  properties  as  the  variation  of  the  flux  density 
with  porosity,  driving  bias,  temperature,  etc.  are  relatively  well  understood.  However,  fluid  flow  through  a 
laboratory  sample  [9,10]  may  differ  from  that  in  an  idealized  percolating  porous  matrix.  Our  goal  is  to 
understand  realistic  porous  media  and  characterize  its  flow'  properties  such  as  transport  of  fluid  and  hydraulic 
conductivity.  The  porous  matrix  is  prepared  in  the  laboratory  from  distribution  of  quartz  sands  with 
appropriate  geomarine  environmental  conditions  such  as  pressure,  compaction,  and  temperature.  Enormous 
efforts  have  been  made  in  recent  years  to  study  structural  patterns  and  flow  of  such  granular  materials  in  a 
variety  of  driving  mechanisms  such  as  vibrations,  mixing,  shearing,  etc.  by  both  computer  simulations  and 
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laboratory  experiments  [11-16].  We  are  focused  on  characterizing  the  flow  through  the  network  of  pores 
formed  from  the  granular  sands,  i.e.,  the  laboratory  sample,  using  a  computer  simulation  model.  Relatively 
few  attempts  [9]  are  made  to  carry  out  such  studies  of  flow  through  such  a  laboratory  porous  medium.  The 
distributions  of  sediments  (barriers)  and  network  of  pores  (pathways)  in  our  laboratory  sample  are  identified 
by  X-ray  computed  tomography  measurements.  Using  a  digitized  image  of  the  sample  as  a  host  matrix, 
attempts  are  made  to  characterize  its  flow  properties  and  hydraulic  conductivity  via  an  interacting  lattice  gas 
method. 

We  consider  a  cubic  lattice  of  size  L  with  L  =  100,  the  digitized  laboratory  sample  with  a  source  of  fluid 
particles  connected  to  the  bottom  of  the  lattice  (z  =  1).  The  sediment  sites  act  as  impenetrable  immobile 
barriers  while  the  nearest-neighbor  pore  sites  attract  the  mobile  hard-core  fluid  particles.  The  interaction 
energy  between  a  particle  and  its  adjacent  pore  site  is  taken  to  be  J  =  —  1.  The  molecular  weight  of  fluid 
constituents  ( M  =  0.1,  relatively  small  here)  facilitates  the  downward  (— z)  fall  and  resist  the  upward  (+z) 
move  via  change  in  the  gravitational  potential  energy  (MAz,  Az  =  ±  1 ).  The  change  in  energy  (interaction  and 
gravitational)  is  used  in  the  Boltzmann  distribution  at  a  fixed  temperature  ( T  =  1)  in  implementing  the  move 
by  the  Metropolis  algorithm.  The  source  of  fluid  particles  at  the  bottom  creates  a  high  concentration  gradient 
which  drives  the  particles  upward  against  gravity.  Additionally,  we  consider  hydrostatic  pressure  bias  (II)  to 
drive  these  constituents  upward  (follows). 

Selection  of  an  adjacent  neighboring  site  along  the  transverse  (±x,  ±y)  and  vertical  (±z)  directions  for  a 
particle’s  move  depends  on  the  bias.  Three  types  of  bias  are  considered  with  0^//^  1  which  co  aid  be  used  to 
describe  different  senario  as  explained  briefly  below. 

I.  Hopping  probability  along  the  transverse  directions  remains  independent  of  bias  which  controls  only  the 
longitudinal  moves.  The  bias  along  the  longitudinal  direction  decreases  with  the  height  (z)  with  a  linear 
gradient. 


P-, 


!-//[!—  (z/L)]  p  _\+H[l-(z/L)] 


(1) 


The  bias  is  linear  with  the  height  z  with  the  largest  at  z  =  0.  Even  at  the  largest  value  of  bias,  H  =  1,  the 
maximum  hopping  probability  along  z-direction  cannot  exceed  j.  Since  the  selection  of  nearest-neighbor  site 
for  hopping  along  the  transverse  directions  (±x,  ±y)  remain  independent  of  bias,  it  is  not  feasible  to 
implement  extreme  bias  that  can  channel  fluid  flow  only  along  one  direction. 

II:  The  hopping  along  the  longitudinal  direction  depends  on  the  bias  but  becomes  independent  of  the  height 
while  maintaining  the  constant  hopping  rate  along  the  transverse  directions. 
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III:  Hopping  along  all  directions  are  coupled  with  bias  such  that  at  extreme  value  of  the  bias  (H  — *  1), 
probability  of  selecting  a  hopping  site  along  the  positive  longitudinal  direction  approaches  one  and  zero  along 
the  rest. 
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With  this  bias,  one  can  examine  a  much  larger  range  in  channeling  the  flow.  Thus,  with  these  three  bias  one 
should  be  able  to  study  flow  in  a  large  pool  of  systems. 

To  implement  the  hopping  [7,8],  we  select  a  particle  say  at  a  lattice  site  i  randomly.  One  of  its  six  nearest- 
neighbor  sites  say  j  is  selected  with  above  probability.  If  site  j  is  empty  then  an  attempt  is  made  to  move  the 
particle  from  site  i  to  site  j  with  the  Metropolis  algorithm.  Periodic  boundary  conditions  are  used  along  the 
transverse  direction;  both  the  top  (z  =  L)  and  the  bottom  (z=  1)  of  the  sample  are  open  for  particles  to 
escape.  A  particle  is  released  from  the  source  into  a  bottom  lattice  site  as  soon  as  it  is  vacated  by  the  particle 
when  it  moves  up  or  diffuses  laterally.  As  usual,  attempts  to  move  each  particle  once  define  the  unit  Monte 
Carlo  step  (MCS)  time.  As  the  simulation  proceeds,  particles  move  in  and  out  of  the  sample.  The  number  of 
particles  and  their  concentration  profiles  (planar  densities  along  x,j>,z)  change  but  eventually  reach  a  steady- 
state — a  non-equilibrium  steady-state  system.  Simulations  are  performed  with  many  independent  initial 
configurations  but  with  the  same  host  matrix,  each  for  a  sufficiently  long  time  steps  to  reach  steady-state. 
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Fig.  1.  Distribution  of  pores  in  the  sample  10(fl  from  quartz  sand  grains  (see  Fig.  2  for  details). 


The  laboratory  sample  is  relatively  dense  with  a  fraction  of  sediment  barriers  pb  =0.718  and  pore  sites 
pQ  =  0.282.  Thus,  it  would  be  better  to  look  at  the  porous  medium  with  empty  sites  of  the  sample.  Fig.  1  shows 
the  distribution  of  pores  in  the  digital  sample.  Note  that  the  fraction  of  pore  sites  is  below  the  site  percolation 
threshold  (pc  0.31 17)  [6],  yet  we  have  connected  pathway  of  pores  that  span  the  entire  sample.  Obviously, 
open  channels  in  the  real  sample  are  different  from  the  idealized  percolating  matrix  with  random,  distribution 
of  barriers  in  a  site  percolation.  To  depict  the  contrast  and  complementarity  in  distribution  of  sediments  and 
pores,  in  Fig.  2  we  present  a  slice  of  the  sample  with  distribution  of  its  sediment  and  pore  sites 

We  perform  the  computer  simulation  experiment  by  injecting  fluid  at  one  end  and  monitoring  its  net  flux 
rate  from  the  opposite  end  by  varying  the  bias  for  three  types  of  driving  mechanisms  described  above.  During 
the  simulation,  we  also  keep  track  of  the  root  mean  square  (rms)  displacement  of  each  particle  (tracer)  and 
their  center  of  mass,  density  profiles,  and  flux  of  fluid.  Fig.  3  shows  the  variation  of  the  rms  displacement  with 
the  time  steps  (?)  for  a  range  of  bias  values  with  pressure  bias  III.  We  see  that  particles,  individually  (tracer) 
and  collectively  (center  of  mass)  drift,  i.e.,  R  oc  t,  in  the  long  time  (asymptotic)  regime.  Such  an  asymptotic 
drift  behavior  is  expected  in  the  long  time  regime  for  our  driven  system  due  to  pressure  bias  and  concentration 
gradient.  Note  that  the  asymptotic  drift  behavior  also  suggests  that  the  system  has  reached  its  steady-state 
which  is  verified  by  the  analysis  of  flux  density  in  the  following. 

There  is  a  net  mass  transfer  due  to  net  driving  field  (along  -Hr  direction)  caused  by  (1)  the  concentration 
gradient  developed  by  the  source  at  the  bottom  and  (2)  the  pressure  bias.  Fluid  particles  enter  at  the  bottom 
and  leave  mostly  from  the  top.  The  amount  of  fluid  entering  the  system  must  be  equal  to  amount  leaving  in  the 
steady-state  due  to  conservation  of  mass,  i.e.,  the  continuity  equation  The  fluid  current  /  can  be  estimated 
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Fig.  2.  Distribution  of  sediment  (left)  and  pores  (right)  in  slices  of  equal  thickness  at  r  =  21-40  in  the  sample  which  is  prepared  from 
quartz  sand  grains.  The  sandy  sediments  are  collected  from  surficial  marine  deposits.  Sediments  are  subrounded  to  subangular  in  shape 
have  a  mean  grain  size  or  d$Q  of  0.375  mm  and  serve  as  the  boundaries  for  the  pores  (right). 
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Fig.  3.  Variation  of  the  root  mean  square  displacement  of  each  particle  (tracer),  R(  and  their  center  of  mass  (right),  Rr  with  bias  III  on 
sample  shown  in  Fig.  1 


from  the  rate  of  net  mass  (Q)  transfer, 
dQ 

/=dT  (4) 

The  current  or  flux  density  (j)  is  the  current  per  unit  cross-sectional  area  ( Ac), 

j  =  i/Ac.  (5) 

The  flux  density  is  linearly  proportional  to  hydraulic  conductivity  (k)  for  flow  with  linear  response  (with 

<5=1), 

j  =  kHs. 


(6) 
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Fig.  4.  Variation  of  the  mass  flux  density  rate  (left)  and  longitudinal  (along  bias)  density  profile  with  bias  III. 
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Fig.  5.  Variation  of  the  steady-state  value  of  flux  density  with  bias  H  for  three  types  1,  II,  Ill. 


Variation  of  the  flux  density  with  the  time  step  is  presented  in  Fig.  4  for  the  entire  range  of  a  bias  (III).  We  see 
that  the  flux  density  reaches  its  constant  value  rather  fast.  Corresponding  density  profiles  are  also  shown  in 
this  figure.  The  non-linear  density  profiles  are  mostly  due  to  morphology  of  the  pore  distribution  except  at 
extreme  bias  H  —  1.  Movement  of  fluid  particles  ceases  as  they  are  trapped  against  sediment  barrier  at  the 
extreme  value  of  type  III  bias. 

Response  of  the  flux  density  to  the  bias  can  be  analyzed  by  examining  the  variation  of  its  steady-state  value 
with  the  bias.  Fig.  5  shows  the  variation  for  three  types  of  pressure  bias.  For  bias  I  and  II,  we  see  relatively 
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linear  responses  of  the  flux  density  to  bias  while  for  bias  III,  a  non-monotonic  dependence  is  cbserved.  Note 
that  in  case  of  bias  I  and  II,  the  hopping  probability  of  fluid  constituents  along  the  transverse  directions  never 
ceases  even  at  extreme  values  of  the  bias.  With  the  pressure  bias  III,  the  movement  in  all  directions  except  r 
ceases  at  H  =  1.  We  believe  that  the  implementation  of  the  pressure  bias  III  is  more  realistic  in  directing  the 
flow  which  shows  a  non-monotonic  response  of  the  flux  density.  We  are  in  a  process  to  carry  out  both 
laboratory  and  computer  simulation  experiments  simultaneously  to  see  which  bias  is  more  appropriate  in 
describing  the  flow  properties  of  the  porous  sediment. 
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